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We study the interaction between two neutral plane-parallel dielectric bodies in the presence of 
a highly asymmetric ionic fluid, containing multivalent as well as monovalent (salt) ions. Image 
charge interactions, due to dielectric discontinuities at the boundaries, as well as effects from ion 
confinement in the slit region between the surfaces are taken fully into account, leading to image- 
generated depletion attraction, ion correlation attraction and steric-like repulsive interactions. We 
investigate these effects by employing a combination of methods including explicit-ion and implicit- 
ion Monte- Carlo simulations, as well as an effective interaction potential analytical theory. The 
latter incorporates strong ion-image charge correlations, which develop in the presence of high 
valency ions in the mixture. We show that the implicit-ion simulations and the proposed analytical 
theory can describe the explicit simulation results on a qualitative level, while excellent quantitative 
agreement can be obtained for sufficiently large monovalent salt concentrations. The resultant 
attractive interaction between the neutral surfaces is shown to be significant, as compared with the 
usual van der Waals interactions between semi-infinite dielectrics, and can thus play a significant 
role at the nano scale. 



I. INTRODUCTION 



Apart from van der Waals (vdW) interactions, no 
universal, and certainly no electrostatic interactions are 
commonly assumed to exist between uncharged (neutral) 
objects [T]-[2. Nevertheless, some interactions are every 
so often observed between neutral surfaces in the pres- 
ence of bathing ionic (salt) solutions [4J and Wernersson 
& Kjellander studied theoretically the case of a salt so- 
lution with mono- and/or divalent ions confined in a slit 
between neutral dielectrics by employing the anisotropic 
hypernetted chain method [5", "F. Electrostatic interac- 
tions, mediated by mobile monovalent ionic solutions, 
are known to be described very well in terms of mean- 
field theories such as the Poisson-Boltzmann (PB) equa- 
tion [71 . For the particular case of symmetric uncharged 
plane-parallel dielectric bodies, which will be considered 
in this paper, the standard PB treatment cannot account 
for inter-surface forces as it yields only a trivial uniform 
ion distribution in the slit between the inner surfaces of 
the slabs, and thus a nil result for the inter-surface force. 
However, in the absence of surface charges, other electro- 
static effects involving mobile ions, not included in the 
PB approach, increase in importance, especially when 
multivalent ions are present in the slit. These include 
fluctuation/correlation and, in particular, dielectric im- 
age effects, which are missing within the standard PB 
description in planar geometry [8]. In fact, correlation 
effects are expected to be significant at uncharged in- 



terfaces, since for each individual multivalent ion, the 
screening cloud and its dielectric images are pinned to 
the mobile ion itself rather than any fixed external sur- 
face charges. Such "ion-image" correlations constitute 
the core insight of our subsequent analysis for neutral sur- 
faces, suggesting that a proper theoretical understanding 
of interactions between neutral surfaces may thus follow 
from the fact that ions and their dielectric images are 
strongly correlated. 

Such correlation effects are expected to be important, 
especially in the presence of highly asymmetric salt mix- 
tures containing a (usually) low concentration of an ionic 
species with a high valency. This latter situation is in fact 
quite relevant experimentally [9,TOJ. We thus focus on 
such asymmetric salt mixtures confined between two un- 
charged plane-parallel dielectric half- spaces, and extend 
the existing analysis [5 to the case where the asymmet- 
ric ionic mixture contains at least one multivalent ion 
component of high valency. The case of symmetric ionic 
fluids and/or weakly charged ionic mixtures at neutral 
interfaces, differing from the focus of the present work, 
has been approached on various levels of approximations 
[3l m [TIHI4] starting from the seminal theoretical work 
of Onsager and Samaras [l5 . 

We base our approach on a combination of methods, 
i.e., explicit-ion, implicit-ion and Yukawa Monte-Carlo 
(MC) simulations as well as an effective "dressed multi- 
valent ion" theory, which accounts for the salt screening 
effects as well as the dielectric polarizability of the neutral 
dielectrics. We show that this latter approach provides 



a simple and useful analytical framework for the present 
and previous ^ results. It is similar to a "dressed counte- 
rion" model discussed in the context of charged surfaces 
in contact with a bathing salt solution containing multi- 
valent "counterions" [16, 17 . There it was shown that 
for high-valency counterions, the degrees of freedom as- 
sociated with the monovalent ions can be integrated out 
and an effective single-particle theory can be obtained 
based on a virial expansion in terms of the counteri- 
ons fugacity, which agrees very well with MC simulations 
within its regime of applicability. This analysis is there- 
fore also akin to the limiting single-particle description, 
known as the strong coupling (SC) theory, which was 
introduced originally for counterion-only systems [18]- 
[25]. The SC theory incorporates, to the leading order, 
strong electrostatic correlations, which develop between 
counterions and the surface charges in the limit of large 
counter ion valency. For uncharged surfaces, such surface 
charge-counterion correlations are absent but, neverthe- 
less, the presence of multivalent ions can make the col- 
lective mean-field description obsolete, provided that the 
surface dielectric ("image charge") effects are properly 
accounted for. As we will show, the single-particle SC 
framework can be reformulated and generalized in order 
to capture the ion-image correlation effects and thereby, 
the interaction mediated by an asymmetric salt mixture 
between two uncharged dielectric bodies, even on a quan- 
titative level. 



These interactions are shown to be quite short-range 
(up to separation distances of just a few nano- meters) 
and are of electrostatic nature even for neutral bodies. 
They combine depletion as well as dielectric image ef- 
fects in a novel kind of interaction that can play a signif- 
icant role in shaping the behavior of neutral soft matter 
systems-most prominently at the nano scale. These in- 
teractions should be therefore considered in parallel with 
the hydrophobic interactions [26] as well as vdW interac- 
tions [2] deemed to play an important role in the inter- 
action of neutral surfaces. 



We base our analysis on explicit-ion MC simulations 
of non-specific ion interactions and while we take dielec- 
tric image interactions properly into account, we disre- 
gard the vdW (dispersion) interactions of the ions with 
the dielectric interfaces, an approximation that is justi- 
fied for high- valency ions [5 . The purpose and empha- 
sis of our endeavor is thus twofold: to show that, for 
highly asymmetric ionic mixtures, a coarse-grained MC 
approach with implicit monovalent ions can be properly 
applied to describe interactions between neutral surfaces, 
and secondly to provide a simple analytical theory, which 
could explain the main features of these simulations on a 
qualitative level that becomes quantitative at sufficiently 
high monovalent salt concentration. 
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FIG. 1: (Color online) Schematic depiction of two neutral 
dielectric half- spaces at separation 2a. The slit region is in 
equilibrium with a bulk reservoir containing a monovalent 1:1 
salt as well as an asymmetric multivalent q:l salt. Explicit- 
ion, implicit-ion and Yukawa levels of description are indi- 
cated schematically (top to bottom)-see the text for details. 



II. MODEL AND METHODS 

Let us consider two neutral plane-parallel dielectric 
half-spaces of dielectric constant e' at a separation dis- 
tance of 2a, as shown schematically in Fig [l] The in- 
tervening slit is assumed to be filled with an aqueous 
solution of dielectric constant £, containing an ionic mix- 
ture comprised of a 1:1 monovalent salt solution of bulk 
concentration no and an asymmetric q:l multivalent salt 
solution of bulk concentration cq, where q stands for the 
charge valency of multivalent ions. The inter facial dielec- 
tric contrast can be quantified by a dielectric mismatch 
parameter^ A = {s — s')/{e -\- e'). In order to investi- 
gate the inter-surface forces induced between neutral di- 
electrics within the model described above, we shall em- 
ploy MC simulations on two different levels of descrip- 
tion as well as an analytical theory. Specifically, we use 
explicit-ion simulations (referred to as "explicit" simu- 
lations) where all ions in the ionic mixture, including 
those of the 1:1 and q:l salt, are simulated explicitly (see 
Fig. [l] top). In this case, the closest-approach distances 
between the ions are taken as = 0.35 nm between 

monovalent cations and anions, and dq^±i = 0.45 nm be- 
tween multivalent ions and monovalent ions. These val- 
ues are small enough so that the excluded- volume effects 
are negligible, yet large enough to prevent Bjerrum pair- 
ing of ions in the solution. The dielectric polarization of 
the half-spaces with planar boundaries is accounted for 
by using repeated image reflections. The image reflection 
sum has to be truncated for computational reasons, and 
we found satisfactory accuracy with only three explicit 
reflections across each surface. All MC simulations were 
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performed using a multi-component grand-canonical en- 
semble, ensuring equilibrium with a bulk solution of the 
prescribed composition. 

In the second approach, we use implicit-ion simula- 
tions (referred to as "implicit" simulations) in which 
only multivalent ions {q > 1) are explicitly simulated, 
whereas the effects of monovalent ions are included im- 
plicitly via a Debye-Hiickel-like ionic cloud. In the im- 
plicit model, apart from the dielectric polarization of the 
half-spaces, we account for the fact that the screening 
ionic cloud around multivalent ions is not spherically 
symmetric but distorted near the (hard) boundaries (see 
Fig. [l] middle). This is accomplished by using the mod- 
ified non-homogeneous Debye-Hiickel (DH) interaction 
kernel (Green's function) iiDH(r, r') in order to compute 
the electrostatic interaction energy between each two (ex- 
plicit) multivalent ions of electric charge qe^ located at 
positions r and in the slit, i.e., (g'eo)^'UDH(r, 1*0- 
Fourier-Bessel representation, ttDH(r, r') can be written 
as 

^dh(Q; z, z') = ii^H(Q; ^0 + ^dh(Q; ^, ^0, (i) 

where Q is the norm of the Fourier wave-vector corre- 
sponding to a two-dimensional Fourier-Bessel transform 
with respect to the lateral directions p = {x^y) perpen- 
dicular to z-axis. The first term above is the isotropic 
(bulk) "Yukawa" potential and the second term, usually 
ascribed to "image" charges, accounts for the presence 
of solid planar dielectric walls. It consistently incorpo- 
rates the dielectric polarization of the half-spaces as well 
as the effects due to confinement of the salt ions to the 
slit only, preventing the spillage of their DH cloud into 
the bounding dielectric half-spaces. This last effect is in 
fact an essential component of the interaction between 
the surfaces. Neglecting this kind of DH cloud distor- 
tions can lead to qualitatively wrong predictions, as we 
shall demonstrate. In explicit form, the Fourier-Bessel 
transforms are given by 
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and k = ^ with k being the Debye screen- 

ing parameter defined via = 47r^B^6- Here £b = 
eQ/{4:7Te£okBT) is the Bjerrum length at ambient tem- 
perature T and 
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is the total density of all monovalent ions in the bulk. 
Note that when multivalent q:l salt is introduced in the 
system, additional monovalent ions with concentration 



qco are simultaneously released into the bulk and thus 
they also contribute to the screening parameter [171 [22] • 
The back-transform of the interaction kernel to real 
space follows from 

1 

UBuir.r') = —- / QdQ Jo{Q\p - p'\)ubu{Q;z,z'), 

^TT Jo 

(6) 

where Jo{x) is the zeroth-order Bessel function of the first 
kind. In the limit of high dielectric mismatch, A ^ 1, the 
interaction kernel can be simplified significantly, yielding 
an analytical form in real space as 
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where = (0, 0, 1) is the unit vector in the +z direc- 
tion and the integer index n runs through both positive 
and negative values. The term with n = stands for the 
direct DH potential, ix^jj(r,r'), of a test charge without 
any dielectric effects. The above expression thus refiects 
a simple picture within which the total potential between 
two parallel dielectric half-spaces can be represented as 
a sum of an infinite number of screened Coulomb image 
charges. This form of image summation significantly re- 
duces also the computational costs, since summing up to, 
e.g., n = 5, already gives a very well-converged result for 
the cases we consider in this work. We have thus used the 
above summation series in the implicit MC simulations, 
and focus on the dielectric jump value A = 0.95, corre- 
sponding to the water-hydrocarbon interface (assuming 
£ = 80 and = 2). The overall Hamiltonian in the 
implicit MC approach then reads, 

H = ^(<7eo)^^DH(ri, r^) + ^ '^{qeo)'^u'^^{Ti, r^), (8) 

i>j i 

combining the screened interaction with distorted DH 
cloud between all multivalent ion pairs and their images 
(first term) as well as their self-image interaction energy 
(second term). 

The implicit MC simulations were performed in the 
grand-canonical ensemble. The simulated lateral sys- 
tem size was chosen such that the slit contained between 
around 50 to 200 multivalent ions. 

In a third simulation approach (referred to as 
"Yukawa" simulations), we use the same implicit model 
as above but neglect the distortion of the ionic clouds 
near the half-spaces as well as their polarizability (see 
Fig. [l] bottom). We thus use the fully symmetric 
screened Coulomb or Yukawa kernel ([2|, which in real 
space reads 
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(9) 



in order to compute the electrostatic interactions between 
multivalent ions. 
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FIG. 2: (Color online) Rescaled density profiles of multivalent ions (q = 4) in the slit between two neutral half-spaces at 
separation 2a — 8nm. The bulk concentration of multivalent 4:1 salt is taken as co = 5 mM and the bulk concentration of the 
additional monovalent salt is no = 30 mM in pane (a) and no = 60 mM in panes (b) and (c). We show the results from our 
explicit MC simulations (orange open squares) , implicit MC simulations (blue filled circles) and the proposed analytical theory 
(solid curves). In (a) and (b) the dielectric mismatch at the boundaries is A = 0.95, whereas in (c) there are no dielectric 
discontinuities, A = 0. The dielectric surfaces are assumed to be impermeable to all ions. 



III. RESULTS 

A. Density profile 

In Fig. [2j we show the density profiles of multivalent 
ions with valency q = 4 within the slit, as obtained from 
our MC simulations for the case of a 4:1 and 1:1 salt mix- 
ture. In all cases the multivalent salt is assumed to have 
a bulk concentration of cq = 5 mM. In Fig. [2^, we have 
an additional 1:1 salt of bulk concentration no = 30 mM 
as well as a dielectric jump parameter of A = 0.95. As 
seen, both explicit and implicit MC simulations show a 
pronounced depletion of multivalent ions from the prox- 
imity of the two neutral walls, where the density tends to 
zero. The explicit and implicit data nearly agree, partic- 
ularly near the walls, with somewhat larger deviations in 
the middle of the slit. A similar surface depletion effect 
can be seen when the additional monovalent salt is in- 
creased but the agreement becomes progressively better; 
in Fig. (2)3, we show the results for no = 60 mM, which 
exhibits a very good agreement between the explicit and 
implicit MC simulations. This confirms the idea that the 
degrees of freedom associated with monovalent ions enter 
only in an effective Debye-Hiickel way as formulated in 
Eqs. Q-Q. Therefore the depletion of multivalent ions 
can indeed be ascribed both to dielectric images due to 
dielectric polarization of the half-spaces as well as to the 
fact that the exclusion of salt ions from the impenetrable 
surfaces leads to a non-trivial distortion of ionic clouds 
in the vicinity of the half-spaces (Fig. [l] middle), effec- 
tively engendering "salt-induced image charges" which 
repel the ions from the walls. Alternatively, one could 
say that ions in the vicinity of the interfaces have fewer 
neighbors and thus all the ions, i.e., both monovalent (not 
shown) [28] and multivalent ions, are depleted from the 
proximity of the surfaces. Such "salt-image effects" are 
present even if there is no dielectric discontinuity at the 
boundaries (A = 0) as seen in Fig.[2]3, and lead to a signif- 



icant depletion near the surfaces, albeit less pronounced 
than in the presence of standard dielectric image effects 
(A > 0, Figs. [2|i and b). Other remarkable differences 
with the latter case are that multivalent ions can have a 
finite contact density at the boundaries when A = and 
also exhibit an even better agreement between explicit 
and implicit MC data (Fig.[2]3). 

The above results clearly go beyond the standard PB 
description which is expected to give only a trivial uni- 
form ion distribution between neutral surfaces in the pla- 
nar geometry. In order to gain any insight into the be- 
havior of this system one needs to conceptualize outside 
of the PB box. We thus proceed by proposing a sim- 
ple analytical theory, akin to the concept of "dressed 
counterions" , which was developed previously in the con- 
text of a multivalent salt solution confined between two 
charged surfaces without a dielectric jump [16] or next 
to one charged dielectric interface [17 , where it was ex- 
haustively and successfully tested against MC simula- 
tions. In the present context with neutral surfaces, where 
one can not speak of "counter ions" proper, we refer to 
this approximation as the "dressed multivalent ion" ap- 
proach. It follows from the same principles as our implicit 
MC simulations, i.e., by integrating out the degrees of 
freedom associated with all monovalent ions treated on 
a linear DH level, and hence by obtaining an effective 
theory where multivalent ions are still described explic- 
itly but the role of monovalent ions enters only implic- 
itly [16 . This effective theory can be justified only for 
highly asymmetric ionic mixtures with large as is of 
primary interest here, and can be obtained systemati- 
cally from a field-theoretical formalism, which we shall 
not discuss further in this paper (see, e.g., Ref. [16 for 
details). Note that the system of dressed multivalent 
ions still presents a full many-body problem and as such 
can be treated exactly only via simulations (which, in 
the present context, would correspond precisely to our 
implicit MC simulations that incorporate the proper im- 



plicit interaction kernels). An approximate analytical 
theory then follows from a virial expansion, which for 
neutral systems is expected to work at sufficiently small 
concentration of multivalent ions, as is often the case in 
experimental systems [9, 10 . For charged surfaces, such 
a virial expansion can be shown to be related to a strong- 
coupling expansion of multivalent counterions [18 , where 
electrostatic correlations between counterions and the op- 
posite surface charges become dominant on the leading 
order. In the present case, such a strong-coupling theory 
for multivalent ions can not be developed as the surfaces 
are neutral and the electrostatic coupling parameter d 
la Netz is non-existent. Nevertheless, the virial expan- 
sion can capture the salient features of the problem as 
it incorporates the ion-image correlation, which in this 
context becomes the dominant mechanism at work for 
multivalent ions at low concentrations. 

On the leading order one thus finds that the only driv- 
ing force for redistribution of multivalent ions in the slit 
\z\ < a is their self-energy, leading to the density profile 
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FIG. 3: (Color online) MC results for the interaction pres- 
sure mediated between two neutral dielectric surfaces as a 
function of their separation in the presence of a multivalent 
4:1 salt mixture of bulk concentration co = 5 mM and with no 
additional monovalent salt (no = 0) and no dielectric discon- 
tinuities (A = 0). We compare the naive Yukawa approach 
with the implicit and explicit MC simulations. 



c(z) = coe-«'™(»'^), 
where the self-energy (per unit /cb^ = reads 



(10) 
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The above analytical result is shown in Fig. ^ (solid 
curves) and appears to capture the main trends of the 
MC data. A close agreement between this theory and 
the simulations is obtained only for large enough mono- 
valent salt concentration (Fig. ^p) or when the dielectric 
jump is small (Fig. In the case A ^ 1, the above 
self-energy expression simplifies into 
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where m = 1,3,5,..., and thus allows for an analogy 
with the Kelvin image charges. 

It should be noted that the above analytical theory pre- 
dicts multivalent ion densities in the slit that are slightly 
smaller than the bulk value cq (see Fig. |2|. In con- 
trast, the simulations do reveal an additional buildup of 
the concentration of multivalent ions in the slit that are 
forced in due to their mutual repulsive interactions in the 
bulk, an effect which is absent in the theory. 



B. Interaction pressure 

We now turn our attention to the interaction (disjoin- 
ing) pressure acting on each of the neutral surfaces due to 
the presence of the salt mixture in the slit. This pressure 
is composed of two separate additive contributions: one 



stemming from the monovalent ions as constituents of the 
1:1 salt as well as of the q:l salt, and the other originating 
from multivalent ions. The latter one can be calculated 
analytically from the single-particle grand canonical po- 
tential within the "dressed multivalent ion" approxima- 
tion. The former is trickier and can only be evaluated 
approximately by starting from an exact pressure equa- 
tion for a simple salt system and then by remorphing it 
in a self-consistent way, using the same line of thought 
as in the case of multivalent ions. This is reasonable as 
the non-homogeneous interaction kernel, Eq. ([T]), applies 
equally to interaction between any two charges, as can 
be clearly deduced from detailed MC simulations [28 . 

In the case of implicit MC simulations as well as the 
proposed analytical approach, the contribution of mono- 
valent ions to the slit pressure, po, can be derived from an 
exact equation formulated by Lovett and Bans [29 , which 
is closely related to the so-called contact-value theorem 
in the plane-parallel geometry [7, 8 , and can be shown 
to have the following approximate linearized (DH) form 

/3po{a) - n(0) - / / dzi dz2 n{zi)n{z2) x 

J_a Jo 

roo Q 

X j Qd(5iXDH(Q;^i,^2) ^^dh(Q;^i,^2)|^_o- (^^) 

It is composed of a van't Hoff term (first term) due to 
the ion osmotic pressure and a correlation term (second 
term) due to the interaction between the ions. It should 
be noted that monovalent ions do also exhibit depletion 
(even though much weaker than the multivalent ions) 
from the vicinity of the walls, and thus a non-uniform 
density profile in the slit. This effect is implicitly em- 
bedded in the modified DH interaction kernel as noted 
above (see also Ref. [13j), and has been verified within 
our explicit MC simulations [28 . We find that the total 
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FIG. 4: (Color online) Pressure acting between the surfaces of two neutral dielectric half-spaces as a function of their separation 
distance for various amounts of salt concentrations, (a) We compare the results from the explicit MC simulations (open squares), 
implicit MC simulations (filled circles), and the proposed analytical theory (solid curves) for a mixture of multivalent 4:1 salt 
with bulk concentration cq — 5mM and different additional 1:1 salt concentrations of no = (black), 30 mM (orange) and 
60 mM (red). In (b) we show the implicit MC results for co = 0.5 mM and various multivalent ion valency q as shown on 
the graph. In (c) we compare the attractive pressures obtained from explicit (squares) and implicit (circles) MC simulations 
with the contribution from the zero-frequency vdW (dashed curves) and full vdW (solid curves) interactions for the cases with 
Co = 5mM and no = 30 mM (black) and 60 mM (red). All our results are on top of the underlying vdW interactions. 



density profile of monovalent ions can be described very 
well by their corresponding self-energy as 



n{z) =n5e-^('''^\ 



(14) 



which is to be used in expression (13) for the depletion 
pressure due to monovalent ions. 

In order to evaluate the slit pressure due to multi- 
valent ions analytically, we now use the virial scheme, 
which to the leading order is given by a single-particle 
grand canonical potential P^'^^^ and can be calculated 
explicitly [T6^, TT . The slit pressure can be deduced by 
differentiation with respect to the volume of the slit, 
Pm = — (S^*^^V^^)/5,A7 where A is the fugacity of the 
multivalent ions. We thus find 
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This contribution should be added to the one from the 
monovalent ions in order to obtain the total slit pres- 
sure. The interaction (disjoining) pressure acting on the 
bounding surfaces follows by subtracting the bulk pres- 
sure, pbuik, as 



p{a) = po{a) + Prnia) - pbu: 



Ik, 



(16) 



where pbuik = Po{cl ^ oo) ^Pm{a ^ oo). 

In order to proceed with the simulation results, we first 
consider the case with no dielectric mismatch across the 
bounding surface, i.e., A = 0. In Fig. |3) we show our 
MC simulation results from the explicit, implicit as well 
as Yukawa approaches (see Section |ll|) for the case where 
the slit contains only a multivalent 4:1 salt mixture of 
bulk concentration cq = 5 mM and no additional mono- 
valent salt (no = 0). As seen, the explicit and implicit 
simulations predict an attractive force at small separa- 
tions which we shall elaborate on further below. It is 



interesting to note however that the naive Yukawa treat- 
ment of the multivalent salt without incorporating the 
image kernel, Eq. ([3|, leads to qualitatively incorrect re- 
sults as it gives a repulsive interaction pressure. The re- 
pulsion in this latter case appears to be purely osmotic in 
nature as surfaces are uncharged. In the case of implicit 
or explicit simulations, the interactions due to dielectric 
and "salt-induced" image charges lead to depletion of 
ions from the slit which in turn engenders a depletion at- 
traction between the surfaces. This depletion is thus gen- 
erated by a confined ionic cloud in contrast to standard 
depletion interactions, which are steric in nature [30] . 
Both types of depletion effects are rather similar except 
that the image-depletion interaction leads to softer de- 
pendence of the interaction pressure on separation when 
compared to the standard hard-sphere depletion. 

In the case where there is a finite dielectric mismatch, 
e.g., A = 0.95, as shown in Fig. |4^, the depletion at- 
traction appears to be more pronounced as more ions are 
depleted from the slit region due both to "salt-induced" 
image as well as dielectric image repulsions. The deple- 
tion attraction increases with the multivalent salt con- 
centration Co (compare Figs. [4^ and b), but also with in- 
creasing the monovalent salt concentration no (Fig. |4^), 
which agrees with the depletion effect seen from the den- 
sity profiles in Fig. [2] At very close separation, all salt 
ions are expelled from the slit and the attractive pressure 
corresponds to the bulk osmotic pressure. In general one 
can discern that the implicit simulations follow nicely the 
explicit simulations, except at very small inter- surface 
separations, where more explicit details on the screening 
cloud is needed. 

The analytical prediction for the total pressure, 
Eq. (16), is shown in Fig. [4^ (solid curves) along with the 
explicit and implicit MC data, demonstrating very good 
agreement with the implicit simulations especially at suf- 
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ficiently large 1:1 added salt concentration. This confirms 
the validity of the single-particle approximation under- 
lying the analytical dressed multivalent ion approach for 
the regime of parameters considered here (i.e., sufficiently 
large sufficiently small multivalent ion concentration 
and sufficiently large monovalent salt concentration). In 
other words, the many-body effects due to interactions 
between different (dressed) multivalent ions is negligible 
and the main contribution results from the self-image 
interactions as included in the analytical theory above. 
The analytical curves show larger disagreement with ex- 
plicit MC data, especially at small separations, but these 
deviations amount to less than about 20% [31 . 

Note that our MC data sometimes reveal a small hump 
at intermediate separations in Fig. [4^. This kind of 
behavior is pronounced at larger valencies and smaller 
screenings, as shown by our implicit MC simulations in 
Fig. |4]3. Comparison with standard steric depletion ef- 
fects [30 leads us to surmise that a larger valency leads 
also to a stronger repulsive interaction between dressed 
multivalent ions that in turn pushes them out of the bulk 
and back into the slit, where they act osmotically, thus in- 
creasing the repulsive component of the interaction pres- 
sure. Note that this "soft sphere packing" mechanism is 
supported by the osmotic pressure measurements on free 
standing films as well as surface force measurements [32j- 
[34j, where an oscillatory behavior is observed between 
neutral (or charged) flat surfaces in the presence of so- 
lutions containing highly charged species (colloidal par- 
ticles or polyelectrolytes) . 

In the final part of this study, we compare the im- 
age depletion interactions discussed above with typical 
values of the vdW contribution that also acts between 
uncharged surfaces in an additive manner, i.e., all the in- 
teractions calculated and analyzed in this work are on top 
of the vdW interactions acting across the slit between two 
dielectric interfaces. In fact the free energy per surface 
area S of the zero- frequency term of the vdW interactions 
(which couples to the charge density fluctuations in the 
salt solution and thus depends on the Debye screening 
parameter) is obtained explicitly as [35j 



f 

Jo 



QdQ 
47r 



ln(l-r2(Q)), (17) 



where T{Q) is defined in Eq. (4|. This term corresponds 
to screened (zero-frequency) electromagnetic field fluctu- 
ations, first addressed by Ninham and Parsegian as de- 
scribed in Refs. [U |3], and later scrutinized at length 
in, e.g., Refs. [5l IH [131 ES]- One should note here 
that in the above free energy we have already taken into 
account a cancellation occurring between the standard 
zero-frequency vdW term and a part of the free energy 
corresponding to an uncoupled, i.e., uncharged, system, 
the result of this cancellation being the above formula. 

At small separations, higher-order Matsubara frequen- 
cies of the electromagnetic field-fluctuations are expected 
to become important as well, whose contribution to the 
interaction free energy of the dielectric half-spaces follows 



from [2] 

s 

with 



n=l 



'QdQ 
2tt 
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[TM TE], 



(18) 



(19) 



where c is the speed of light in vacuo, and e(z^n) and 
IJi{iin) are the dielectric response function and the mag- 
netic permeability of the slit at imaginary Matsubara fre- 
quencies, respectively. TM and TE correspond to trans- 
verse magnetic and transverse electric modes and the 
n summation is over imaginary Matsubara frequencies 
= 27rnkBT/h^ where h is the Planck constant divided 
by 27r. For the sake of simplicity, we assume that for 
all layers fi{t^n) — 1 and that the dielectric function for 
the slit is that for water. ATM(^^n) and ATE(^^^n) are 
functions of the frequency-dependent response functions 
of the semi-infinite dielectrics and the water slab whose 
explicit forms can be found in Ref. [2j. The non-zero 
Matsubara frequency vdW free energy can be written as 



>1 



487ra2 ' 



(20) 



where A^>i is the corresponding Hamaker coefficient. 
We take An>\ = 3 zJ, which is on the upper bound for 
the non-zero Matsubara frequency contribution for hy- 
drocarbons interacting across an aqueous medium [37] . 
The complete Hamaker coefficient A^=o + A^>i for, e.g., 
dimyristoyl phosphatidylcholine (DMPC) and dipalmi- 
toyl phosphatidylcholine (DPPC) multilayers has a typ- 
ical value in the range of 2.87 — 9.19 zJ [38 . 

As seen in Fig. [4]3, the pressure due to the image de- 
pletion interactions is of the same order of magnitude as 
the total vdW (but larger than the zero-frequency vdW 
contribution) at separations distances about 1 nm and 
above, especially at large monovalent salt concentrations. 
One should not forget, however, that the vdW interac- 
tions and the effects discussed in this work are additive, 
hence the total interaction pressure is the sum of the 
vdW pressure and Eq. (16). Obviously the image-related 



interactions discussed here play as significant a role at 
the nano scale as the vdW interactions and should thus 
not be overlooked. 



IV. CONCLUSION 

We consider an asymmetric ionic fluid, consisting of 
multi-species monovalent salt ions and a single species of 
multivalent ions, confined between two uncharged (neu- 
tral) dielectric half-spaces. Specifically we studied the 
ion distribution and the electrostatic interaction medi- 
ated by the ionic fluid between two bounding dielectric 
interfaces. 
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Apart from vdW interactions, the standard DLVO the- 
ory of colloidal interactions in an ionic solution [2 [71 18] 
predicts that electrostatic interactions between neutral 
dielectrics are non-existent. Any interaction found be- 
sides the vdW interactions would thus present a seri- 
ous challenge to accepted wisdom in the colloidal domain 
(such electrostatic forces are reported between randomly 
charged dielectrics that are neutral on the average but 
carry a charge disorder component j39H43] , which we do 
not consider here and focus only on strictly uncharged 
dielectrics). 

Our purpose was twofold: first, to investigate the na- 
ture of the interactions engendered by the dielectric in- 
homogeneities and salt redistribution between the dielec- 
tric interfaces by means of both explicit- and implicit-ion 
simulations, and to assess the validity of a simple analyt- 
ical description that allows us to couple in a consistent 
way both the "salt-induced" image effects as well as the 
dielectric polarization effects. We showed that the ionic 
mixture of a monovalent salt and a salt with multivalent 
ions can be decomposed into two parts, which can be 
treated in a single analytical framework. We described 
the monovalent salt solution on a DH level, which-in 
other words-means that the monovalent ions are inte- 
grated out, leading to a simplified but accurate descrip- 
tion for the multivalent ions in the solution. We refer to 
this procedure as "dressed multivalent ion" approxima- 
tion. This approximation is valid for highly asymmetric 
salt mixtures, i.e., when multivalent ions have a large 
charge valency [16l [17] . Even with the monovalent ion 
degrees of freedom integrated out, the remaining system 
containing only the dressed multivalent ions still repre- 
sents a full many-body problem and additional approxi- 
mations need to be introduced in order to calculate the 
partition function of the system analytically. To this end, 
we employed a virial expansion in order to obtain the dis- 
tribution of multivalent ions and the interaction pressure 
between the surfaces mediated by them up to the leading 
order in multivalent counterion fugacity. This amounts 



to a simple single-particle theory for dressed multivalent 
ions which is expected to work for sufficiently small mul- 
tivalent ion concentration and high enough monovalent 
ion concentration as we have indeed confirmed by our 
explicit-ion simulations. 

Although these approaches are based on different ratio- 
nales, the final outcomes share some qualitative similar- 
ities. All the ions are repelled from neutral surfaces due 
to "salt-induced" image effects, which are always present 
for ion-impermeable interfaces, as well as dielectric im- 
age effects present only in the case of dielectric inhomo- 
geneities. This leads to ionic depletion in the slit, which 
results in an attraction between the slabs caused by the 
osmotic pressure of multivalent ions in the bulk. Our the- 
oretical results are corroborated by explicit and implicit 
MC simulations, which show an excellent agreement in 
the cases where a comparison makes sense. 

The magnitude of the effects described in this work 
is comparable or surpasses the underlying vdW inter- 
actions between neutral dielectric interfaces, and should 
thus perforce be included in colloidal interaction equi- 
libria for uncharged matter in complex multicomponent 
ionic solutions. 
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